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Motivated by recent experiments on volborthite single crystals showing a wide |-magnetization 
plateau, we perform microscopic modeling by means of density functional theory (DFT) with the 
single-crystal structural data as a starting point. Using DFT-|-t/, we find four leading magnetic ex¬ 
changes: antiferromagnetic J and J 2 , as well as ferromagnetic J' and Ji. Simulations of the derived 
spin Hamiltonian show good agreement with the experimental low-field magnetic susceptibility and 
high-held magnetization data. The |-plateau phase pertains to polarized magnetic trimers formed 
by strong J bonds. An effective J —>■ 00 model shows a tendency towards condensation of magnon 
bound states preceding the plateau phase. 

PACS numbers: 71.70.Gm, 75.10.Jm, 75.30.Et, 75.60.Ej 


The perplexing connection between quantum mag¬ 
netism and topological states of matter renewed inter¬ 
est in frustrated spin systems [1]. A prime example is 
the S' = 5 antiferromagnetic kagome Heisenberg model 
(KHM), whose ground state (GS) can be a gapped topo¬ 
logical spin liquid, as suggested by large-scale density- 
matrix renormalization group (DMRG) simulations [2, 3]. 
Although DMRG results were recently corroborated by 
nuclear magnetic resonance (NMR) measurements on 
herbertsmithite [4] , alternative methods vouch for a gap¬ 
less spin liquid [5, 6] and the discussion is still not settled. 

One of the remarkable properties of the KHM is the 
presence of field-induced gapped phases that manifest 
themselves as magnetization plateaus [7-9] . A key ingre¬ 
dient thereof are closed hexagonal loops of the kagome 
lattice that underlie the formation of valence-bond solid 
states [8]. By far widest is the ^-magnetization plateau, 
whose structure is well described by singlets residing on 
closed hexagons, and polarized spins (Fig. 1, left) [7-10]. 

Despite the considerable progress in understanding 
both quantum and topological aspects of the KHM, most 
theoretical findings still await their experimental verifi¬ 
cation. The reason is the scarceness of material real¬ 
izations: only a handful of candidate KHM materials is 
known to date. A prominent example is herbertsmithite, 
where S = ^ spins localized on Cu^'*" form a regular 
kagome lattice [11]. Other candidate materials feature 
exchange couplings beyond KHM as kapellasite [12-17], 
haydeeite [12, 16-19], francisite [20], or barlowite [21, 22]. 

The natural mineral volborthite Gu3V207(0H)2-2H20 
was considered a promising KHM material [24, 25] , until 
it was noticed that the local environment of two crys- 
tallographically distinct Gu sites hints at different mag¬ 
netically active orbitals [26]. Density functional theory 



FIG. 1. (Color online) The structure of the |-magnetization 
plateau in the kagome model (KHM), coupled frustrated 
chains (CFG) model from Ref. [23], and the J-J'-Ji-J 2 model. 


(DFT) calculations show that this has dramatic impli¬ 
cations for the spin physics, giving rise to coupled frus¬ 
trated chains (CFG) with ferromagnetic (FM) nearest- 
neighbor and antiferromagnetic (AF) second-neighbor 
exchanges, and interstitial spins that are AF coupled 
to the two neighboring chains [23]. However, detailed 
structural studies reveal that below ~300 K all Cu atoms 
have the dx 2 _y 2 as the magnetically active orbital [27], 
questioning the applicability of the CFG model for vol¬ 
borthite. Furthermore, the CFG model features the 
I-magnetization plateau with a semiclassical “up-up- 
down” structure (Fig. 1, middle), which was never ob¬ 
served in powder samples [28, 29]. Recent magnetization 
measurements on single crystals overturned the experi¬ 
mental situation; a broad ^-magnetization plateau sets 
in at i7ci —26T and continues up to at least 74 T [30]. 

Puzzled by the remarkable difference between the 
single-crystal and powder data, we adopt the structural 
model from Ref. [30] and perform DFT and DFT-l-17 cal¬ 
culations. We find a microscopic model which is even 
more involved than CFG: besides sizable Ji and J2 form¬ 
ing frustrated spin chains, the coupling between the chain 
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FIG. 2. (Color online) (a) Microscopic magnetic model of 
volborthite featuring four relevant exchange couplings: an¬ 
tiferromagnetic J (thick bars) and J 2 (solid curved lines), 
as well as ferromagnetic J' (dashed lines) and Ji (wiggly 
lines). Magnetic trimers formed by J exchanges are high¬ 
lighted (shaded ovals). Magnetic Cu atoms are shown as 
large spheres within CUO 4 squares, nonmagnetic V atoms are 
middle-sized spheres within VO 4 tetrahedra. (b) The Cu- 
O-V-O-Cu superexchange paths in the magnetic trimer. (c) 
Magnetic trimers form a basis for (d) the effective model with 
ferromagnetic J \, as well as antiferromagnetic J 2 , J 2 , and Jz. 

and the interstitial Cu atoms is now facilitated by two in¬ 
equivalent exchanges, a sizable AF J and a much weaker 
FM J' . Due to the dominance of J, the magnetic planes 
break up into magnetic trimers (Fig. 2). By using exact 
diagonalization (ED) of the spin Hamiltonian, we demon¬ 
strate that this model agrees with the experimental mag¬ 
netization data and explains the nature of the plateau 
phase (Fig. 1, right). Further insight into the low-field 
and low-temperature properties of volborthite is provided 
by analysing effective models of pseudospin-^ moments 
T living on trimers. Thus, a model based on effective ex¬ 
changes Ji, J 2 and supports the presence of a bond 
nematic phase due to condensation of two-magnon bound 
states. Finally, we conjecture that powder samples of 
volborthite suffer from disorder effects pertaining to the 
stretching distortion of Cu octahedra. 

We start our analysis with a careful consideration of 
the crystal structure. Volborthite features a layered 
structure, with kagome-like planes that are well sepa¬ 
rated by water molecules and non-magnetic V 2 O 7 groups. 
Magnetic Cu^'*' atoms within the planes occupy two dif¬ 
ferent sites: Cu(2) with four short Cu-0 bonds forms 
edge-sharing chains, and interstitial Cu(l) located in be¬ 
tween the chains. Different structural models in the 


TABLE I. Direct Cu..Cu distances dcu..Cu (in A), transfer in¬ 
tegrals t (in meV) and exchange integrals J (in K). GGA-|-[/ 
results are provided for three different values of the on-site 
Goulomb repulsion Ud- The two numbers in each entry per¬ 
tain to the two structurally inequivalent layers; this minor 
layer dependence is ignored in the subsequent analysis. 

. , J (GGA-ht/) 

^Cu..Cu ^ 

_ Ud = 8.5 eV 9.5 eV 10.5 eV 

J 3.053/3.058 -191/-194 193/205 156/167 127/136 

J' 3.016/3.020 -80/-84 -29/-22 -30/-25 -32/-26 

Ji 2.922/2.923 -98/-100 -65/-65 -76/-74 -77/-76 

J 2 5.842/5.842 64/64 32/31 26/22 22/21 


literature suggest either squeezed [31] or stretched [32] 
Cu(l)06 octahedra. The DFT study of Ref. [23] em¬ 
ployed a structure with a squeezed Cu(l) octahedron. 
Although such configuration can be realized at high tem¬ 
peratures [27], Cu(l)06 octahedra are actually stretched 
in the temperature range relevant to magnetism [27, 30]. 
The respective structural model was never studied with 
DFT, hence we fill this gap with the present study. 

For DFT calculations [33], we use the generalized gra¬ 
dient approximation (GGA) [34] as implemented in the 
full-potential code fplo9. 07-41 [35]. We start with a 
critical examination of all structural models proposed so 
far, by optimizing the H coordinates and comparing the 
total energies. In this way, we find that the single crystal 
structure of Ref. [30] has the lowest total energy [33] . All 
further calculations are done for this structural data set. 

To evaluate the magnetic couplings, we project the 
relevant GGA bands onto Cu-centered Wannier func¬ 
tions [33]. The leading transfer integrals t (>50meV) of 
the resulting one-orbital (dx^-y^) model are provided in 
Table I. Their squared values are proportional to the AF 
superexchange, which is usually the leading contribution 
to the magnetism. However, such one-orbital model fully 
neglects FM contributions that are particularly strong for 
short-range couplings (dcu..Cu ^ 3 A). Hence, to evaluate 
the exchange integrals that comprise AF and FM contri¬ 
butions, we perform DFT-I-C/ calculations for magnetic 
supercells and map the total energies onto a Heisenberg 
model. These results are summarized in Table I. 

Prior to discussing the magnetic model, we should note 
that the structural model of Ref. [30] implies the pres¬ 
ence of two similar, albeit symmetrically inequivalent 
magnetic layers, with slightly different Cu..Cu distances. 
Since the respective transfer (t) and exchange (J) inte¬ 
grals for both layers are nearly identical (Table I), we can 
approximately assume that all layers are equal and halve 
the number of independent terms in the model. 

The resulting four exchanges, J, J', Ji, and J 2 form 
the 2D microscopic magnetic model depicted in Fig. 2. 
This model is topologically equivalent to the GFG model: 
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it consists of chains with first- (Ji) and second-neighbor 
(J 2 ) couplings and the interstitial Cu atoms coupled to 
two neighboring chains. However, the exchange between 
the interstitial spins and the chains is realized by two 
different terms: a dominant AF J and much weaker FM 

. This contrasts with the CFG model, where both ex¬ 
changes are equivalent {J = J'). 

From the structural considerations, the difference be¬ 
tween J and J' may seem bewildering, as Cu..Cu dis¬ 
tances (Table I) and Cu-O-Cu angles (104.6° versus 
102.4°) are very similar. Indeed, for the usual Cu-0- 
Cu path, the superexchange would be only marginally 
different for J and J'. The difference originates from 
the long-range Cu-O-V-O-Cu path (Fig. 2, b) which 
provides an additional contribution to J, but not J', 
since the latter lacks a bridging VO 4 tetrahedron. It 
is known that long-range superexchange involving empty 
V d states can facilitate sizable magnetic exchange of up 
to 300 K [36]. Hence, it is the long-range Cu-O-V-O-Cu 
superexchange that renders J much stronger than J^ 

A distinct hierarchy of the exchanges J > | Jij > J 2 , J' 
leads to a simple and instructive physical picture. The 
dominant exchange J couples spins into trimers that tile 
the magnetic layers. Each trimer is connected to its four 
nearest neighbors by FM J' and Ji , and to its two second- 
neighbors by AF J 2 (Fig. 2). In contrast to the CFG 
model, where frustration is driven exclusively by J 2 , the 
coupled trimer model has an additional source of frustra¬ 
tion: triangular loops formed by J, J' and Ji. Together 
with J 2 , they act against long-range magnetic ordering. 

DFT-|-C/-based numerical estimates for the leading ex¬ 
change couplings allow us to address the experimen¬ 
tal data. To simulate the temperature dependence of 
the magnetic susceptibility Xj ED of the spin Hamil¬ 
tonian is performed on lattices of N = 2A spins, us¬ 
ing the approximate ratios of the exchange integrals 
J:J':Ji:J 2 = 1:—0.2:—0.5:0.2 (Table I). The simulated 
curves are fitted to the experiment by treating the overall 
energy scale J, the Lande factor g and the temperature- 
independent contribution xo as free parameters. In this 
way, we obtain a good fit down to 35 K with J = 252 K, 
(7 = 2.151 and xo = l-06x 10“'^ emu / [mol Gu] (Fig. 3). 
ED even reproduces the broad maximum at 18 K, which 
stems from short-range antiferromagnetic correlations. 
Deviations at lower temperatures are finite-size effects. 

After establishing good agreement with the x(^) data, 
we employ a larger lattice of N = 36 spins and calcu¬ 
late the GS magnetization curve, which shows a wide 
|-magnetization plateau between the critical fields iJci 
and Hc 2 (Fig. 3, bottom left). Scaling with J and 
g from the xi^) fit, without any adjustable parame¬ 
ters, yields Hd = 22 T in agreement with the experimen¬ 
tal i?ci=26T. In the plateau phase, first- and second- 
neighbor spin correlations within each trimer amount to 
(So • Si) = (Si • S 2 ) = -0.4938 and (Sq • S 2 ) = 0.2470, very 
close to the isolated trimer result (—^ and respec- 



FIG. 3. (Color online) Top: magnetic susceptibility of the 
microscopic spin Hamiltonian calculated by ED on a A = 24 
site lattice compared to experiment (Ref. [37]) and an iso¬ 
lated trimer model. Bottom left: GS magnetization curve 
simulated on a lattice of A=36 spins for the same model. In¬ 
sets are magnifications of the respective data. Bottom right: 
GS magnetization of the full effective model [33] with N = 24, 
26, and 30 pseudospins compared to experiment (Ref. [30]). 

tively [33]). Hence, the ^-plateau phase can be approxi¬ 
mated by a product of polarized spin trimers formed by 
strong J bonds (Fig. 1, right), and thus is very different 
from the plateau phases of the KHM (Fig. 1, left) and 
the GFG model (Fig. 1, middle). The plateau stretches 
up to a remarkably high — 225 T, at which the spin 
trimers break up, allowing the magnetization to triple. 

ED-simulated spin correlations indicate that the sim¬ 
plest effective model — the isolated trimer model — 
already captures the nature of this plateau phase. On 
general grounds, we can expect the isolated trimer model 
to be valid only at high temperatures. However, it pro¬ 
vides a surprisingly good fit for magnetic susceptibility 
down to 60K (Fig. 3, top), i.e. at a much weaker energy 
scale than the leading exchange J ~250K (Fig. 3). This 
motivates us to treat the inter-trimer couplings perturba- 
tively and derive a more elaborate effective model valid 
at low temperatures and in low fields (T, gg^H/k^ J). 

To this end, we adopt the lowest-energy doublet of each 
trimer at 77 = 0 as the basis for a pseudospin- ^ operator 
Ti. The i-plateau phase corresponds to the full polariza¬ 
tion of pseudospin-i moments (m®® = ™sat/3). Degen¬ 
erate perturbation theory to second order in the inter- 
trimer couplings yields an effective Heisenberg model on 
a triangular lattice with spatially anisotropic nearest- 
neighbor couplings Ji = —34.9K and J2=36.5K and 
much weaker longer-range couplings such as 1/2 = 6.8 K 
and = 4.6 K shown in Fig. 2 (d) [33]. The competition 
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between FM Ji and AF underlies the frustrated na¬ 
ture of the effective model. Larger finite lattices available 
to ED of the effective model allow us to amend the crit¬ 
ical field iJci estimate compared to the full microscopic 
model (Fig. 3, bottom left) and reproduce a pronounced 
change in the M{H) slope (Fig. 3, bottom right), which 
agrees with the experimental kink at ~22T [30]. 

Effective models provide important insights into the 
nature of field-induced states. Recent NMR experiments 
on single crystals revealed the emergence of the incom¬ 
mensurate collinear spin-density-wave (SDW) phase “IF 
(iJ<23T) and the “N” phase preceding the plateau 
(23 T < iL < 26 T) [30]. We first address the nature of the 
latter phase, by treating the fully-polarized pseudospin 
state (the 1/3-plateau state of volborthite) as the vac¬ 
uum and analyzing the magnon instabilities to it. 

To this end, we resort to a model with three lead¬ 
ing effective couplings J7i, Ji, and J!^- This model is 
equivalent to the frustrated FM square lattice model, 
where a bond nematic order emerges owing to con¬ 
densation of two-magnon bound states (bimagnons) for 
Ji = J 2 > 0.4] J7i I [38]. Here we take the approximate ra¬ 
tio ^ 72 / 1 ^ 1 ] = 1 of the perturbative estimates, and study 
the influence of J 2 on the ground state. We find that 
the bond nematic order is robust for J 2 l\Ji \ 0.3 [33], 

as signaled by the occurrence of bimagnon condensation 
at iL/ib at which the plateau state is already destabi¬ 
lized, but before single-magnon condensation sets in at 
(Fig. 4, a). The bond nematic phase shows no long- 
range magnetic order besides the field-induced moment, 
but it is characterized by a bond order with an alternat¬ 
ing sign of directors Dij = {T^T^ — T^Tj) residing on Ji 
bonds (Fig. 4, b) [39]. This phase is a viable candidate 
for the experimentally observed “N” phase, whose NMR 
spectra are not explained by simple magnetic orders [30] . 

While bimagnons are stable in a wide region of the Ji- 
J1-J2 model, longer-range effective couplings such as J 3 
tend to destabilize bimagnons. However, a slight tun¬ 
ing of the microscopic model (e.g., increasing of |J'|) 
can counteract this effect, thereby recovering the nematic 
phase [33]. Long-range effective couplings are also sen¬ 
sitive to weak long-range exchanges neglected in the full 
microscopic model. In the absence of experimental esti¬ 
mates for these small exchanges, the J\-Ji-J 2 effective 
model is an adequate approximation, which allows us to 
study the nature of the field-induced phases in volbor¬ 
thite. 

Below 23 T, NMR spectra indicate the onset of an 
incommensurate collinear phase “H” [30]. Unfortu¬ 
nately, incommensurate spin correlations produce irregu¬ 
lar finite-size effects that impede an ED simulation. Yet, 
on a qualitative level, further truncation of the model to 
the effective couplings J\ and Ji leads to an anisotropic 
triangular model, for which a field-theory analysis pre¬ 
dicts the SDW order for m < = |msat [40]. 



FIG. 4. (Color online) (a) The behavior of one- and two- 
magnon gaps in the |-plateau phase, which gives rise to a 
bond nematic phase, (b) Schematic picture of the bond ne¬ 
matic phase in the effective model. Orientation of dark el¬ 
lipses represents the sign of directors Dij on J\ bonds [39]. 


Next, we go a step beyond the Heisenberg model 
and consider antisymmetric Dzyaloshinskii-Moriya (DM) 
components for the leading couplings J and Ji. 
By performing noncollinear DFT-I-C/ calculations with 
vasp [41], we obtain jHij/Ji ~0.12 with Di nearly or¬ 
thogonal to the frustrated chains. DM vectors D within 
the trimers are nearly orthogonal to the respective in¬ 
teratomic vectors and amount to |Dj/J~0.09 [33]. We 
analyzed the influence of D for isolated dimers with ED 
and found a minute change in spin correlations in the 
plateau state, which amounts to 2% at most. However, 
these DM interactions are the leading anisotropy at low 
fields, and can give rise to the two consecutive transitions 
to the incommensurate phase “I” (T < 1K, iJ < 4T) [27]. 

Finally, we address the intriguing question why the | 
plateau has not been observed in powder samples. We 
remind the reader that the trimers are underlain by the 
stretching distortion of Cu(l )06 octahedra, which selects 
two out of four neighboring VO 4 octahedra for the J 
superexchange pathway (Fig. 2, b). In single crystals, the 
distortion axes are fixed, and the trimers form an ordered 
parquet-like pattern. Powder samples on the other hand 
are more prone to a random choice of the distortion axis. 
A single defect of this type permutes J and J', ruining the 
trimer picture locally. This tentative scenario explains 
the absence of a plateau and the strong dependence on 
the sample quality in the powder magnetization data. 

In summary, the stretching distortion of the magnetic 
Cu(l )06 octahedra in volborthite leads to the model 
of coupled trimers, very different from the anisotropic 
kagome and coupled frustrated chain models discussed 
in earlier studies. Based on DFT calculations and ED 
simulations, we conclude that i) the microscopic mag¬ 
netic model of volborthite contains four exchanges with a 
ratio J: J': Ji:J 2 = 1:—0.2:—0.5:0.2 and J = 252K; ii) the 
I-magnetization plateau can be understood as a prod¬ 
uct of nearly independent polarized trimers, and iii) the 
effective J1-J1-J2 model shows indications for a bond 
nematic phase which precedes the onset of the plateau. 

Note added: A recent NMR study [42] supports our 
bond nematic phase scenario below the |-plateau. 
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Supplemental Material for Magnetic Behavior of Volborthite Cu 3 V 207 ( 0 H) 2 - 2 H 20 
Determined by Conpled Trimers Rather than Prnstrated Chains 

O. Janson, S. Furukawa, T. Momoi, P. Sindzingre, J. Richter, and K. Held 

TABLE SI. Comparison of total energies yielded by nonmagnetic DPT calculations using the full-potential code fplo version 
9.07-41 [1]. The respective references, experimental method (XRD for x-ray diffraction, ND for neutron diffraction), number of 
formula units in the unit cell (Z), the space group and the fc-mesh are given in columns 1-5. and correspond to the 

total energies obtained using parameterizations for the exchange and correlation potential from Refs. [2] and [3], respectively. 
The hydrogen positions are optimized with respect to total energy, the resulting forces are below 0.01 eV/A. 


reference 

method 

Z 

sp. gr. 

fc-mesh 

EGGAjz (Hartree) 

E^GKjz (Hartree) 

Ref. [4] 

XRD 

1 

C2/m 

8x8x8 

-7695.23828 

-7678.79575 

Ref. [5] 

ND 

1 

C2/m 

8x8x8 

-7695.24387 

-7678.80053 

Ref. [5] 

XRD 

1 

C2lm 

8x8x8 

-7695.22987 

-7678.79266 

Ref. [6] 

XRD 

2 

CC 

4x4x4 

-7695.24642 

-7678.80498 

Ref. [7] 

XRD, 150 K 

2 

C2/c 

4x4x4 

-7695.24236 

-7678.80240 

Ref. [7] 

XRD, 322 K 

1 

C2lm 

8x8x8 

-7695.23947 

-7678.79349 

Ref. [8] 

XRD 

2 

C2/c 

4x4x4 

-7695.24878 

-7678.80835 

Ref. [9] 

XRD 

4 

P2i/c 

4x4x4 

-7695.25191 

-7678.81252 


TABLE S2. Crystal structure used for DPT and DPT-l-17 calculations. We use the 50 K data from Ref. [9] in the conventional 
setup, i.e. the space group is P2i/c (14), a =14.41 A, 6 = 5.8415 A, c = 10.6489 A, and /3 = 95.586°. The hydrogen coordinates 
shown bold are optimized within GGA-fP. All other atomic coordinates are taken from Ref. [9]. 


atom 

Wyckoff position 

xja 

J//6 

zjc 

Cul 

2 d 

1 

2 

0 

1 

2 

Cu2 

2 c 

0 

0 

1 

2 

Cu3 

4e 

0.00233 

0.75600 

0.24626 

Cu4 

4e 

0.49930 

0.74625 

0.25280 

VI 

4e 

0.12660 

0.52334 

0.49667 

V2 

4e 

0.37297 

-0.02131 

0.00304 

oi 

4e 

0.24991 

-0.04890 

0.00158 

02 

4e 

0.05938 

0.00418 

0.34368 

HI 

4e 

0.12961 

0.01212 

0.34767 

03 

4e 

0.55895 

-0.00105 

0.34372 

H2 

4e 

0.62959 

0.00012 

0.34817 

04 

4e 

0.07440 

0.50571 

0.34159 

05 

4e 

0.42524 

-0.00532 

0.15831 

06 

4e 

0.26060 

0.50190 

0.17836 

07 

4e 

0.24083 

0.03230 

0.31923 

08 

4e 

0.09912 

0.20580 

0.07522 

09 

4e 

0.60273 

0.78700 

0.07584 

OlO 

4e 

0.08846 

0.73814 

0.07118 

on 

4e 

0.41442 

0.75383 

0.42937 

H3 

4e 

0.20595 

0.56499 

0.12695 

H4 

4e 

0.23942 

0.35162 

0.20575 

H5 

4e 

0.26676 

0.88701 

0.29281 

H6 

4e 

0.29312 

0.10203 

0.37325 












TABLE S3. Comparison of transfer integrals tij (meV) evaluated using Wannier functions for one-orbital (a:^ — , only) and 

two-orbital {x^ — y^ and 3z^ — r^) models. For the latter, only the transfer integrals between the half-filled — y^ are provided. 
The respective band dispersions are shown in Fig. 3 and Fig. SI. Direct Cu..Cu distances dcu..Cu are given in A). 


transfer integral 

dcu..Cu 

t for one-orbital WFs 

t for two-orbital WFs {x^ — y'^ x^ — y'^) 

t 

3.053/3.058 

-191/-194 

-198/-195 

t' 

3.016/3.020 

-80/-84 

-80/-85 

ti 

2.922/2.923 

-98/-100 

-98/-100 

t2 

5.842/5.842 

64/64 

66/66 



o GGA 

- WF 2-orbitals (c/, 2 .y 2 and model 



rx SY rzu rt z 

/f-vector 


FIG. SI. (Color online) Left: GGA band structure of volborthite (solid lines) with band characters for the Cu d^ 2 _y 2 orbitals. 
The dashed line indicates the projection onto Cu-centered Wannier functions (WFs) within a one-orbital x^ — y^ model. Note the 
deviations around 0.4 eV. Right: The deviations around 0.4eV are remedied in a two-orbital {x^ —y^ and 3z^ — r^) model. 
The fc-points are: r=(000). X= (^OO), M= (^fO), Y= (OfO), Z= (OO^), XZ= (fOf), MZ= (f f f), and YZ= (Of f). 


EFFECTIVE HAMILTONIANS 

In the main text, we analyze the physics at low temper¬ 
atures and low fields (T, h = gfisH/kB ^ J) hy using an 
effective Hamiltonian of pseudospin- ^ moments living on 
trimers. Here we describe details of the derivation of this 
effective Hamiltonian and our analyses of the magnon 
spectra and the magnetization process of this model. We 
also derive yet another low-temperature effective Hamil¬ 
tonian suitable for high fields (T, |h — 3 J/2| <C J), which 
can be used to determine the high-held end Hc 2 of the 
plateau phase and the saturation held iJc 3 - 

Our use of these effective Hamiltonians is motivated 
by the distinct hierarchy of the exchange couplings J > 
|Ji| > J 2 , |d'| revealed in DFT-(-17. We start from the 
limit of large J, where the system decouples into trimers. 
In this limit, each isolated trimer exhibits a wide ^- 
magnetization plateau over 0 < h < 3J/2. At each 
endpoint of the plateau, where a level crossing occurs 
in the ground state of each trimer, the total system pos¬ 
sesses a macroscopic degeneracy and is highly susceptible 
to inter-trimer couplings. Our effective Hamiltonian is 
derived around each endpoint by performing degenerate 


perturbation theory in terms of J', Ji, and J 2 . This kind 
of approach is known as strong coupling expansion [10], 
and is often used to analyze quantum spin systems with 
coupled cluster structures such as coupled dimers [10- 
12] and coupled trimers [13, 14]. Our derivation of the 
effective Hamiltonians goes essentially in parallel with 
Refs. [13, 14]. 


Isolated trimer 

We first consider an isolated trimer Hamiltonian 

2 

Hoi2 = J{o-Si + S^-S2)-hYySl ( 1 ) 

i=o 

where Si and S' 0,2 are the spin-^ operators at the central 
and other sites of the trimer, respectively. At h = 0, 
because of the SU(2) symmetry and the parity symmetry 
around the site 1, the eigenstates are classified into the 
quadruplet {Jg^i)}, the even-parity doublet and 

the odd-parity doublet {IdJ^)}, where fi is the eigenvalue 
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FIG. S2. The spectrum of an isolated trimer Hamiltonian 
(1) as a function of the magnetic field h. The eigenstates and 
eigenenergies are shown in Eqs. (2) and (3), respectively. 


of Their wave functions for /r > 0 are given by 

|g+3) = I ttt), 

k+i) ^ (I 

M+i) ^ (I 

lTi> = ;^(im)- 

Other eigenstates with /r < 0 are obtained by apply¬ 
ing the spin reversal to the above. The Zeeman term in 
Eq. (1) commutes with the Heisenberg terms, and only 
shifts the eigenenergies by —hfi. The eigenenergies in the 
presence of a field h are calculated as 

|9m) : \ j 

\d^) - J - h^i■ ( 3 ) 

\d'f,) ■■ - hfi, 

which are plotted in Fig. S2. At h = 0, the ground states 
are the even-parity doublet |(i±i); they are split for h ^ 
0. At h = 3J/2, a level crossing occurs, and the ground 
state is replaced by the fully polarized state l^+a). This 
level structure gives a wide magnetization plateau with 
'm/msat = 1/3 over 0<h<3J/2. At the level crossing 
points h = 0 and 3J/2, the total system consisting of Nt 
trimers possesses a macroscopic degeneracy 2^‘ of ground 
states. We restrict ourselves to this degenerate manifold, 
and analyze the splitting of the degeneracy due to inter- 
trimer couplings by deriving effective Hamiltonians. 

The small Hilbert space of the isolated trimer model 
facilitates a direct evaluation of its partition function, 
and hence its magnetic susceptibility, which is given by: 


I itt) +1 nt)), 

I nt) - 2 | nt)), 
int)). 



FIG. S3. Effective Heisenberg Hamiltonain for |/i| <C J. The 
five largest couplings obtained in the second-order perturba¬ 
tion theory are displayed. For each trimer indicated by a light 
red bar (with its center at r) , we introduce a pseudospin-1 op¬ 
erator Tr- These pseudospins form a triangular lattice. This 
lattice consists of two sublattices A and B corresponding to 
two different directions of trimers. The vectors u and v (of 
lengths |u|, |n| ~ b ~ 6A) connect between neighboring sites 
on the triangular lattice. 


X*(/3)=/3(^ + ^(2 + exp^ +exp^) (4) 

where j5 is the inverse temperature. 


Effective Hamiltonian for |m/msat| <1/3 

Here we derive the effective Hamiltonian for the range 
|w/msat| < 1/3. We label each trimer by the position r of 
its central site. These positions form a triangular lattice 
as shown in Fig. S3. This lattice consists of two sublat¬ 
tices A and B corresponding to two different directions 
of trimers. 

For \h\ <C J, we use the lowest-energy doublet |d-i-i)r 
at h = 0 as the local basis on each trimer r. Using these 
states, we introduce a pseudospin-1 operator 



where a = are the Pauli matrices. 

The first-order effective Hamiltonian is derived by pro¬ 
jecting the inter-trimer interactions onto the degenerate 
manifold Vq — Span({|dj_i )r}). The obtained Hamil¬ 
tonian is a Heisenberg model of pseudospins {Tr} on 
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TABLE S4. Nonzero coupling constants in the effective Hamiltonian (7). The vectors u and v are defined in Fig. S3. 

Only for J-i and the relative vector Ar depends on the sublattice X, as indicated in the corresponding rows. The first- 
and second-order perturbative estimates (third and fourth columns) are calculated by substituting Eq. (9) into Eqs (6) and 
(8). We also consider Models I and II shown in the fifth and sixth columns, where the leading three and four couplings in 
the second-order model are taken into account, respectively; these simplified models help to understand the essential physics 
arising from the leading couplings. For each of these models, one- and two-magnon condensation points, /ks 

with n = 1, 2, are presented. Here, g = 2.151 obtained in the main text is used to convert into In the four models, 

the one-magnon spectrum has the minimum at incommensurate wave vectors kb = ±(<3,0), and the value of QI{2ti) is also 
presented. 



relative vectors Ar 

Ist-order 

2nd-order 

Model I 

Model H 

Jl 

U, V 

-44.8 K 

-34.9 K 

-34.9 K 

-34.9 K 

J 2 

U + V 

44.8 K 

36.5 K 

36.5 K 

36.5 K 

Ji 

— U -b V 


6.8 K 

6.8 K 

6.8 K 

Ji 

{2u,A), {2v,B) 


4.6 K 


4.6 K 

Ji 

{2v,A), {2u,B) 


1.7 K 



Ji 

2 u + v,u + 2u 


1.7 K 



Js 

2 (u + v) 


-1.3 K 



f.;;* {hS) 


22.4 K (15.5 T) 

35.0 K (24.2 T) 

19.9 K (13.8 T) 

27.2 K (18.9 T) 

Q/{2-k) 


0.333 

0.370 

0.341 

0.360 

(HiV) 


22.4 K (15.5 T) 

35.0 K (24.2 T) 

25.7 K (17.8 T) 

29.0 K (20.1 T) 


the triangular lattice with spatially anisotropic nearest- 
neighbor couplings 


^(Ist) _ 



n 





( 6 ) 


as shown in Fig. S3. 

To calculate the second-order effective Hamiltonian, 
we take into account various processes where a state 
in Vo is virtually promoted to an excited state out¬ 
side Vq, and then comes back to Vq by the operations 
of the inter-trimer interactions. Such virtual processes 
give rise to various further-neighbor interactions between 
pseudospins. The resulting Hamiltonian is a Heisenberg 
model 


HeS = EE *J^Ar,Xr'^r ‘ ^r+Ar ^ E^- (7) 

r Ar r 


order expressions are given by 


^^(2nd) 


^( 2 nd) 

J 2 


„/(2nd) 

J 2 


^( 2 nd) 

'-'3 


^/(2nd) 

'-'3 


^( 2 nd) 

J 4 


^( 2 nd) 

'-'5 


r/^ 211Ji2 -b48JiJ' - 118J'^ 

= ) +-ifflOj- 

8J2(-4Ji -b 5J') 

243J ’ 

8 J2^ 2(-2Ji -b J')(-13Ji -b 8J') 

~ ~ 243J 

2(-2Ji -b -8J') 

“ 243J ’ 

2(-5Ji -b4J')(-Ti -4J') 

“ 243J ’ 

5(-2 Ji -b J')^ 

486J 

8J2(-4Ji -b 5J') 

243J 

32J2" 

243J ■ 

( 8 ) 


Most of these interactions do not depend on the sub¬ 
lattice X, and have the translational invariance of the 
triangular lattice. Only J 3 and depend on the sub¬ 
lattice X, and double the unit cell of the effective model 
when Jo ^ J^', the primitive vectors of the system then 
change to ±m -b v. Using 


J : J' : Ji : J 2 = 1 : -0.2 : -0.5 : 0.2, J ~ 252 K, (9) 


where X^ = A, B indicates the sublattice which r belongs 
to. The coupling constants J\r,x show seven nonzero 
different values as listed in Table S4, and their second- 


obtained in the main text, the effective coupling con¬ 
stants are calculated as in the columns “Ist-order” and 
“2nd-order” in Table S4. The five largest couplings in 
the second-order model are displayed in Fig. S3. 
















11 


Magnon spectra 


state forms a single band 


Using the effective model TJeff in Eq. (7), we here cal¬ 
culate the spectra of one- and two-magnon excitations 
which lower the total magnetization by one and two, re¬ 
spectively, from the |-magnetization plateau state. The 
plateau state corresponds to the fully polarized state of 
pseudospins M+i)r =: |vac), which we view as the 
magnon vacuum in the following. Then, := ± iT^ 

play the roles of magnon annihilation and creation op¬ 
erators, and the magnon occupation number at the site 
r is given by := 5 — Tjf. Using these operators, the 
Hamiltonian is rewritten as 

HeS = Nt + h'^^Tlr + Hmag, ( 10 ) 

where 

r Ar 

= 2j7i + 1/2 + + 2^/4 + iTs, (11) 

^mag ^ ^ 

r 

r Ar 

( 12 ) 


e{k) = -J -P J\r cos(fc ■ Ar). (13) 

Ar 

When Ja ^ the unit cell is doubled, and the one- 
magnon bands are obtained by diagonalizing the 2x2 
matrix 

M{k) = eo{k)I + J(k) ■ a, (14) 

where I is the identity matrix, and 

eo(fc) = - J + Ji cos{ku + ky) + j!^ cos(A:„ - ky) 

+ + vi^a) [cos(2fc„) -|- cos(2A:^,)] 

+ A cos{2ky -\- 2ky), 

J^{k) = Ji {cosky-\-cosky) (15) 

“t” [cos(2/i;.,i -t- ky) -\- cos{ky -t- 2A;^j)], 

jy{k) = 0 , 

= ^(*^3 - J3) [cos(2fc„) - cos(2fc«)] 

with ky = k ■ u and ky = k ■ v. The two bands are 
calculated as 

e±(fc) = eo(fc) ± I J(/c)|. (16) 


2(j^r+?^r+Ar+h.C.) TlyTly-^/^y 


The first term in Eq. (10) is the energy of the vacuum 
I vac). The second term is the Zeeman term, which plays 
the role of a magnon chemical potential. The third term 
Hmag contains magnon kinetic and interaction terms, and 
does not depend on h. Below we calculate the n-magnon 
spectrum of i7mag, and in particular determine the lowest 
energy in it. We then find from Eq. (10) that the 
minimum energy cost for creating n magnons from the 
vacuum is given by + hn. This energy reaches zero at 
h = /n, which corresponds to an n-magnon 

condensation point. 

While the second-order effective model is expected to 
describe well the properties of the original model with 
Eq. (9) for T,h J, it contains small couplings of mag¬ 
nitudes less than 2 K, which may be influenced easily by 
potential parameter changes or inclusion of small further- 
neighbor couplings in the original model, which are not 
estimated in the DET-I-U calculation. We therefore con¬ 
sider also Models 1 and II shown in Table S4, where 
the leading three and four couplings in the second-order 
model are taken into account, respectively. These sim¬ 
plified models help to understand the essential physics 
arising from the leading effective couplings. 

We first analyze one-magnon spectra of i7mag- When 
= J 3 , the effective model has the translational invari¬ 
ance of the triangular lattice. In this case, one-magnon 


The (lower) one-magnon band e{k) or e_(/c) is plot¬ 
ted for the four models in Fig. S4. In all the cases, the 
spectrum has the minimum energy E^^'^ at incommensu¬ 
rate wave vectors kb = ±(Q,0). The one-magnon con¬ 
densation point and the incommensurate 

wave number QI{2Tr) are presented in Table S4. Remark¬ 
ably, the band e{k) for Model I shows a line of nearly 
degenerate minima extending roughly along the vertical 
direction. This implies suppression of magnon hopping 
(and hence relative enhancement of magnon interactions) 
along this direction. 

We next analyze two-magnon spectra. To this end, we 
introduce the two-magnon basis [15], which is given by 

|Ar, A:) = ^ ^ e^'^-TUT.-^Jvac) (17) 
when ^73 = ^ 7 ' 3 , and by 

1 ^’ = - 7 ^ T. (18) 

with X = A,B when jTa 77 Here, the magnon relative 
vector Ar is chosen in the range |Ar |/6 < O(IO^) in such 
a way that the double counting of ±Ar which lead to the 
same state is avoided. The matrix elements of Hmag are 
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(a) 1 st-order 


180 


(b) 2nd-order 



-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k)(b / 271 

(c) Model I 


140 


C\J 


>. 



-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k^b / 271 


0.6 

0.4 


0.2 


- 0.2 

-0.4 


- 0.6 



i-i 40 


- 20 


- -20 


-40 


-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k^b / 271 

(d) Model II 


0.6 


0.4 


0.2 


- 0.2 


-0.4 


- 0.6 



- 40 


- 20 


- 0 


- -20 


-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k^b / 271 


FIG. S4. One-magnon spectra of Jlmag in Eq. (12) for the four models shown in Table S4. We plot e{k) in Eq. (13) for (a) and 
(c) , and e-{k) in Eq. (16) for (b) and (d). The color indicates an energy in units of Kelvin. A hexagon or a square indicates 
the first Brillouin zone. The wave numbers, kx and ky, are defined along the horizontal and vertical directions of Fig. S3, 
respectively. 


given by 

(Ar, fc|iJi„ag|Ar',fc) 

= ^ e*'=-’'(vac|r+T+,i/„agT;-T,-^,,|vac), (19) 

r 

{X,Ar,k\H^,^\X',Ar',k) 

^ik-(r-uSxB) 

r^X' 

X (vac|Tj^^^Tj^^^_i_^j,iJi„agJ(. T(.+at-'(20) 

for the bases in Eqs. (17) and (18), respectively. We note 
that the dependence on Nt has dropped in these expres¬ 
sions, and we are effectively treating an infinite system 
(an error can only arise from the finite cutoff for |Ar|). 
By performing Lanczos diagonalization for such matri¬ 
ces, we have calculated the lowest eigenenergy for each k, 
which is plotted in Fig. S5. In (a) the first-order and (b) 
second-order models, we find nearly degenerate minima 
at kb = (0,0) and ±(2(5,0) with energy 
these can be interpreted as two independent magnons. 


each with the one-magnon lowest energy In (c) 

Model I and (d) Model II, by contrast, a single minimum 
at kb = (0, 0) with energy E^‘^'> < 2E^^'> is formed, imply¬ 
ing a significant effect of attractive interactions. In (c) 
Model I, in particular, two magnons acquire an apprecia¬ 
ble binding energy 2E^^'> — E^^'^ = 11.6 K, which is likely 
to be due to aforementioned suppression of magnon hop¬ 
ping along the vertical direction. The two-magnon con¬ 
densation point = —E^‘^^j2 is presented in Table S4. 
The relation seen in Models I and II indicates 

that bimagnon condensation leading to a bond nematic 
order occurs with lowering the field h from the plateau 
phase. We have also calculated 3- and 4-magnon spec¬ 
tra (not shown), finding no indication of multimagnon 
condensation with h!'^^ > h!'^^ {n = 3,4). 
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(a) 1 st-order 



-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k)(b / 231 

(c) Model I 



-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k)(b / 231 


(b) 2nd-order 



(d) Model II 



-0.6 -0.4 -0.2 0 0.2 0.4 0.6 

k^b / 231 


FIG. S5. Two-magnon spectra of flmag in Eq. (12) for the four models shown in Table S4. 


Exact diagonalization 

In the main text, we present the exact diagonalization 
result for the magnetization process of the second-order 
effective model in Eq. (7). Here we explain the cluster 
shapes used for this study. A finite-size cluster is specified 
by two vectors Lj = L^jU + LyjV [j = 1, 2; L^j, (L^j -I- 
Lyj)j2 G Z], which set the periodic boundary conditions 
Tr = Tr+Lj ■ The number of pseudospins in the system 
is given by Nt = We have used the 

following clusters specified by {Lui,Lyi) and (T« 2 ,T„ 2 ): 

iVt = 24: (5,3), (-3,3); 

iVt = 26: (5, !),(-!, 5); (21) 

A^t = 30: (5,-1), (0,6). 


choices, making the cluster shape as isotropic as possi¬ 
ble. For example, 120° rotation of (Lui,L„i) = (5,1) 
gives (—— Lyi) = (—1,4). However, the latter 
vector connects between different sublattices, and thus 
we shift it slightly and set {Lu 2 , Ly 2 ) = (—Ij 5) to obtain 
the Nj = 26 cluster. Exact diagonalization of the effec¬ 
tive model was performed using TITPACK ver. 2 [16]. 

Effective Hamiltonian for 1/3 < m/maat < 1 

For |/i —3J/2| «C J, we can derive yet another effective 
model valid in the range 1/3 < m/rusat < 1, by using 
\q^3)r and |(i+i)r as the local basis on each trimer r. 
Using these states, we introduce a pseudospin-^ operator 


The reason for these choices is as follows. In finite-size 
clusters with periodic boundary conditions, two distinct 
interactions Ty ■ T^+Ar and Ty ■ Ty^/j^yi fall into the iden¬ 
tical one if Ar — Ar' = miLi + TO 2 T 2 ('^ 11,3712 G ^)- 
In the above choices, where Li and L 2 are sufficiently 
long, such a situation can be avoided. Furthermore, 
Li and L 2 form nearly 120° to each other in the above 


fy = 



f r(9-|-§ l\ 


( 22 ) 


The all down state of the pseudospins (T* = —1/2) cor¬ 
responds to the |-magnetization plateau of the original 
model; the all up state (T/ = -1-1/2) corresponds to the 
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saturation. The first-order effective Hamiltonian is de¬ 
rived in a way similar to the previous case, and has the 
form of an XXZ model on the triangular lattice with spa¬ 
tially anisotropic couplings. The Hamiltonian is given by 




r Ar *- 


^xy/r^xrj^x \ 

-^r-\-Ar -^r “^r+Ar/ 


I rrz rpzrpz _ 1 rpZ 

^ c/Ar-^r-^r+Ar '^ / . ^ 


(23) 


where 


j-y := j-v = j-y = -iJi- 2J'), 

_ 1 

lL-\-V 2 

*^2 •— J^u-\-v — T7^'^2 


_ ^xy _ Z 7 

^2 ' u^v o ^^ 2 5 


(24) 


18 

h = h — ho, Hq = —J y^(5Ti -|- 5 J 2 -l- IIT^). 

Z lo 


Using Eq. (9), the parameters in this model are calculated 
as 


= (- 4 . 2 ,- 9 . 1 ,16.8, 2 . 8 ) K, 

ho = 326 K. 


Because of the spin-reversal symmetry of the XXZ Hamil¬ 
tonian, the magnetization process of the original model 
is symmetric about h = hp at this order of perturbation 
theory. 

We determine the saturation field hg of the effective 
model (23) by analyzing the single-magnon instability of 
the saturated state. One-magnon excitation band above 
the saturated state is calculated as 

e(fc) =h - 2J^ - J 2 + Jf*'(cos fc„ -I- cos K) 

+ J2^ COs(fc„ -I- fc„). 

When 0 < < 23^^, this has the minimum at kb = 

I I 

±((3,0) with Q = 2arccos and the saturation field 

hs is determined as the held h at which this minimum 
reaches zero, hence 

/ ■rxy\2 

hs = ± ± . (27) 


Using Eq. (25), the high-held end h^i of the plateau 
and the saturation held /ics of the original model are 
calculated as 


(/ic2, /Jcs) = (^0 - K-, ho + hs) = (324,328) K 
[(i?c2,ifc3) = (224,227) T]. 

These helds agree well with the held at which rapid in¬ 
crease of mjmsat is found in the exact diagonalization 
result of the original model. This rapid increase in the 
short held width Hco — Hc 2 = 2hs = 2.6 T arises from the 
rather small coupling constants [Eq. (25)] in the effective 
XXZ model (23). For Hc 2 < H < Hcs, the system is 
expected to exhibit an incommensurate magnetic order 
with a wave vector Q/{2'k) = 0.46. 


* olegjanson@gmail.com 
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